# Load libraries and data ----------------------------------------------------------
library(tidyverse)
library(estimatr)
library(texreg)
library(nnet)
library(AER)
library(kableExtra)
library(fixest)
library(multcomp)
library(lmtest)

source("make_table.R")
df1 <- read_csv("experiment1_cleaned.csv")
df2 <- read_csv("experiment2_cleaned.csv")
dfA <- read_csv("experimentA_cleaned.csv")
dfB <- read_csv("experimentB_cleaned.csv")
dfC <- read_csv("experimentC_cleaned.csv")
dfC_respondents <- read_csv("experimentCresp_cleaned.csv")

# set path, which is location for output, here
# for example:
# save_path <- "C:/PublicReplicationFile/output"

options(knitr.table.format = "latex")
select <- dplyr::select

# * Experiment 1 ----------------------------------------------------------

# Table A1. Sample demographics -----------------------------------------------

exp1_sample_df <- df1 %>% 
  select(O_general08, O_general10, O_general11, O_general12,
           female, gender_miss, african_am, white, minority_unknown, age_imputed, age_miss)

exp1_sample <- as.data.frame(round(apply(exp1_sample_df, 2, mean, na.rm = TRUE), digits = 3)) %>% 
  cbind(round(apply(exp1_sample_df, 2, sd, na.rm = TRUE), digits = 3)) %>% 
  cbind(round(apply(exp1_sample_df, 2, min, na.rm = TRUE), digits = 3)) %>% 
  cbind(round(apply(exp1_sample_df, 2, max, na.rm = TRUE), digits = 3)) %>% 
  rownames_to_column() %>% 
  mutate(rowname = case_when(rowname == "O_general08" ~ "Voted in 2008",
                             rowname == "O_general10" ~ "Voted in 2010",
                             rowname == "O_general11" ~ "Voted in 2011",
                             rowname == "O_general12" ~ "Voted in 2012",
                             rowname == "female" ~ "Female",
                             rowname == "gender_miss" ~ "Missing Gender",
                             rowname == "african_am" ~ "Black",
                             rowname == "white" ~ "White",
                             rowname == "minority_unknown" ~ "Missing Race",
                             rowname == "age_imputed" ~ "Age (imputed)",
                             rowname == "age_miss" ~ "Missing Age"))

make_kable(exp1_sample, col.names = c("Variable", "Mean", "SD", "Min", "Max"),
           caption = "Summary Statistics, Experiment 1 (MS)")

# Table A2. Balance testing ----------------------------------------------------------------
ml1 <- multinom(treatment ~ O_general12 + O_general11 + O_general10 + O_general08 + 
         female + gender_miss + african_am + white + minority_unknown + age_imputed + age_miss, data = df1)

null_ml1 <- multinom(treatment ~ 1, data = df1)

lrtest(ml1, null_ml1)

ml1_extracted <- texreg::extract(ml1, beside = TRUE,
                         include.aic = FALSE, include.bic = FALSE,
                         include.loglik = FALSE, include.deviance = FALSE,
                         include.groups = FALSE)

make_texreg("exp1_balancetests", list(ml1_extracted[[1]], ml1_extracted[[2]],
                                      ml1_extracted[[3]], ml1_extracted[[4]]),
            custom.model.names = c("Gloating Villain", "Foiled Villain",
                                   "Happy Hero", "Disappointed Hero"),
            custom.coef.map = list("O_general08" = "Voted in 2008",
                                   "O_general10" = "Voted in 2010",
                                   "O_general11" = "Voted in 2011",
                                   "O_general12" = "Voted in 2012",
                                   "female" = "Female",
                                   "gender_miss" = "Missing Gender",
                                   "african_am" = "Black",
                                   "white" = "White",
                                   "minority_unknown" = "Missing Race",
                                   "age_imputed" = "Age (imputed)",
                                   "age_miss" = "Missing Age",
                                   "(Intercept)" = "Constant"),
            custom.note = paste("\\item[\\hspace{-5mm}] %stars.",
                                "\\item[\\hspace{-5mm}] \\textit{Note:} Models estimated using multinomial logistic regression.
                                The dependent variable is each treatment condition regressed on the demographic covariates 
                                used in the main analysis of Experiment 1. The p-value of the likelihood ratio test 
                                was 0.94, which means we fail to find that the covariates 
                                jointly predict treatment assignment."),
            multiple.tasks = "glm",
            scalebox = .8,
            caption = "Balance Tests of Treatment Conditions, Experiment 1 (MS)")
  




# Table 3. Experiment 1 -----------------------------------------------------------------
lm1 <- lm_robust(voted14_gen ~ treat_gloatingvillain + treat_foiledvillain + treat_happyhero + treat_disappointedhero,
                data = df1)

lm2 <- lm_robust(voted14_gen ~ treat_gloatingvillain + treat_foiledvillain + treat_happyhero + treat_disappointedhero + 
                   O_general12 + O_general11 + O_general10 + O_general08 + 
                   female + gender_miss + white + african_am + minority_unknown + age_imputed + 
                   age_miss,
                 data = df1)


make_texreg("exp1_regression", list(lm1, lm2),
            custom.model.names = c("Turnout in 2014", "Turnout in 2014"),
            custom.coef.map = list("treat_gloatingvillain" = "Gloating Villain",
                                   "treat_foiledvillain" = "Foiled Villain",
                                   "treat_happyhero" = "Happy Hero",
                                   "treat_disappointedhero" = "Disappointed Hero",
                                   "O_general08" = "Voted in 2008",
                                   "O_general10" = "Voted in 2010",
                                   "O_general11" = "Voted in 2011",
                                   "O_general12" = "Voted in 2012",
                                   "female" = "Female",
                                   "gender_miss" = "Missing Gender",
                                   "african_am" = "Black",
                                   "white" = "White",
                                   "minority_unknown" = "Missing Race",
                                   "age_imputed" = "Age (imputed)",
                                   "age_miss" = "Missing Age",
                                   "(Intercept)" = "Constant"),
            custom.note = paste("\\item[\\hspace{-5mm}] %stars.",
                                "\\item[\\hspace{-5mm}] \\textit{Note:} Models estimated using ordinary least squares regression, with robust standard errors. 
                                The dependent variable is a binary measure of turnout in 
                                the 2014 Mississippi Special Election.
                                Baseline of untreated Control group (versus the 4 experimental treatment groups), Male (versus Female and Gender Unknown), and
                                White (versus Black and Race Unknown/Other). Along with the dummy variable ``Age Unknown'' indicating when the age covariate is missing values,
                                we replace the missing values in ``Age (imputed)'' with the mean age."),
            multiple.tasks = "normal",
            scalebox = .8,
            caption = "Effect of Gloating Villain on Turnout, Experiment 1 (MS)")




# Analysis for Experiment 1 -----------------------------------------------
lm1$coefficients[1]

lm1$coefficients[2]/lm1$coefficients[1]

lm2$coefficients[4]
lm2$coefficients[5]
lm2$coefficients[3]

summary(glht(lm2, linfct = c("treat_gloatingvillain = 0")))
summary(glht(lm2, linfct = c("treat_gloatingvillain - treat_disappointedhero = 0")))
summary(glht(lm2, linfct = c("treat_gloatingvillain - treat_foiledvillain = 0")))
summary(glht(lm2, linfct = c("treat_gloatingvillain - treat_happyhero = 0")))
summary(glht(lm2, linfct = c("treat_foiledvillain - treat_happyhero = 0")))
summary(glht(lm2, linfct = c("treat_foiledvillain - treat_disappointedhero = 0")))
summary(glht(lm2, linfct = c("treat_happyhero - treat_disappointedhero = 0")))

# *Experiment 2 ------------------------------------------------------------

# Table 4. Experiment 2 -----------------------------------------------------------------
lm1_base_field2 <- feols(voted ~ treat_villain + treat_reportcard, cluster = ~vpc_hhid,
                         data = df2)

lm1_all_field2 <- feols(voted ~ treat_villain + treat_reportcard + racenum + female + 
                          age + age2 + married + voted2012g + voted2014g + voted2016g + voted2018g + 
                          bc_hh_size + statehousedistrict + catalistmodel_ideology + ideology2, 
                        cluster = ~vpc_hhid,
                        data = df2)

lm1_h7_field2 <- feols(voted ~ treat_villain + treat_reportcard + racenum + female + 
                         age + age2 + married + voted2012g + voted2014g + voted2016g + voted2018g + 
                         bc_hh_size + catalistmodel_ideology + ideology2, 
                       cluster = ~vpc_hhid, subset=~statehousedistrict == 7,
                       data = df2)

lm1_h38_field2 <- feols(voted ~ treat_villain + treat_reportcard + racenum + female + 
                          age + age2 + married + voted2012g + voted2014g + voted2016g + voted2018g + 
                          bc_hh_size + catalistmodel_ideology + ideology2, 
                        cluster = ~vpc_hhid, subset=~statehousedistrict == 38,
                        data = df2)
make_texreg("exp2_regression", list(lm1_base_field2, lm1_all_field2, lm1_h7_field2, lm1_h38_field2),
            custom.model.names = c("All ", "All", "\\makecell{Florida House\\\\District 7}", 
                                   "\\makecell{Florida House\\\\District 38}"),
            custom.coef.map = list("treat_villain" = "Gloating Villain",
                                   "treat_reportcard" = "Report Card",
                                   "racenumasian" = "Asian",
                                   "racenumblack" = "Black",
                                   "racenumhispanic" = "Hispanic",
                                   "racenumother" = "Race Other/Unknown",
                                   "female" = "Female",
                                   "age" = "Age",
                                   "age2" = "Age^2",
                                   "married" = "Married",
                                   "voted2012g" = "Voted in 2012",
                                   "voted2014g" = "Voted in 2014",
                                   "voted2016g" = "Voted in 2016",
                                   "voted2018g" = "Voted in 2018",
                                   "bc_hh_size" = "Household Size",
                                   "catalistmodel_ideology" = "Catalist Ideology",
                                   "ideology2" = "Catalist Ideology^2",
                                   "statehousedistrict38"= "Florida House District 38",
                                   
                                   "(Intercept)" = "Constant"),
            custom.note = paste("\\item[\\hspace{-5mm}] %stars.",
                                "\\item[\\hspace{-5mm}] \\textit{Note:} Models estimated using ordinary least squares regression, 
                                with standard errors clustered by household. The dependent variable is a binary measure of turnout in 
                                the 2019 Florida Special Election.
                                Baseline of untreated Control group (versus the Gloating Villain and Report Card treatment groups),
                                Male (versus Female and Gender Unknown),
                                White (versus Asian, Black, Hispanic, and Other), 
                                Male/Gender Unknown (versus Female), and House District 7 (versus House District 38)."),
            multiple.tasks = "feols",
            reorder.gof = c(2,1),
            scalebox = .8,
            caption = "Effect of Gloating Villain on Turnout, Experiment 2 (FL)")



# Analysis for Experiment 2 -----------------------------------------------
lm1_all_field2$coefficients[2]

lm1_base_field2$coefficients[2]/lm1_base_field2$coefficients[1]

lm1_all_field2$coefficients[3]

lm1_h7_field2$coefficients[2]-lm1_h38_field2$coefficients[2]

summary(glht(lm1_all_field2, linfct = c("treat_villain - treat_reportcard = 0")))

# Interaction between state house district and each treatment

interaction <- feols(voted ~ treat_villain*statehousedistrict + treat_reportcard*statehousedistrict + racenum + female + 
                       age + age2 + married + voted2012g + voted2014g + voted2016g + voted2018g + 
                       bc_hh_size + catalistmodel_ideology + ideology2, 
                     cluster = ~vpc_hhid,
                     data = df2)
summary(interaction)


# Table B1. Sample demographics -------------------------------------------

exp2_sample_df <- df2 %>% 
  select(voted2012g, voted2014g, voted2016g, voted2018g, white,
         asian, black, hispanic, female, age, married, bc_hh_size,
         catalistmodel_ideology) |> 
  mutate(catalistmodel_ideology_white = case_when(white == 1 ~ catalistmodel_ideology,
                                                  TRUE ~ NA_real_)) |> 
  mutate(catalistmodel_ideology_minority = case_when(white == 0 ~ catalistmodel_ideology,
                                                     TRUE ~ NA_real_)) |> 
  select(-catalistmodel_ideology)

exp2_sample <- as.data.frame(round(apply(exp2_sample_df, 2, mean, na.rm = TRUE), digits = 3)) %>% 
  cbind(round(apply(exp2_sample_df, 2, sd, na.rm = TRUE), digits = 3)) %>% 
  cbind(round(apply(exp2_sample_df, 2, min, na.rm = TRUE), digits = 3)) %>% 
  cbind(round(apply(exp2_sample_df, 2, max, na.rm = TRUE), digits = 3)) %>% 
  rownames_to_column() %>% 
  mutate(rowname = case_when(rowname == "voted2012g" ~ "Voted in 2012",
                             rowname == "voted2014g" ~ "Voted in 2014",
                             rowname == "voted2016g" ~ "Voted in 2016",
                             rowname == "voted2018g" ~ "Voted in 2018",
                             rowname == "female" ~ "Female",
                             rowname == "black" ~ "Black",
                             rowname == "white" ~ "White",
                             rowname == "asian" ~ "Asian",
                             rowname == "hispanic" ~ "Hispanic",
                             rowname == "age" ~ "Age",
                             rowname == "married" ~ "Married",
                             rowname == "bc_hh_size" ~ "Household Size",
                             rowname == "catalistmodel_ideology_white" ~ "Catalist Ideology (White)",
                             rowname == "catalistmodel_ideology_minority" ~ "Catalist Ideology (Non-White)"))

make_kable(exp2_sample, col.names = c("Variable", "Mean", "SD", "Min", "Max"),
           caption = "Summary Statistics, Experiment 2 (FL)")
# Table B2. Balance testing ----------------------------------------------------------------
## Create multinomial logistic regression for FL HD 7
ml2_7 <- multinom(condition ~ bc_elections_voted + bc_mailaddressuspstype_H + age + 
                black + hispanic + asian + female + married + bc_hh_size,
                data = df2, subset = statehousedistrict == 7)

## Calculate null model, then calculate chi-squared statistic and p-value
null_7 <- multinom(condition ~ 1, data = df2, subset = statehousedistrict == 7)

lrtest(ml2_7, null_7)


## Create multinomial logistic regression for FL HD 7
ml2_38 <- multinom(condition ~ bc_elections_voted + bc_mailaddressuspstype_H + age + 
                    black + hispanic + asian + female + married + bc_hh_size,
                  data = df2, subset = statehousedistrict == 38)

## Calculate null model, then calculate chi-squared statistic and p-value
null_38 <- multinom(condition ~ 1, data = df2, subset = statehousedistrict == 38)

lrtest(ml2_38, null_38)

## Make regression tables for 
ml2_7_extracted <- texreg::extract(ml2_7, beside = TRUE,
                         include.aic = FALSE, include.bic = FALSE,
                         include.loglik = FALSE, include.deviance = FALSE,
                         include.groups = FALSE)

ml2_38_extracted <- texreg::extract(ml2_38, beside = TRUE,
                           include.aic = FALSE, include.bic = FALSE,
                           include.loglik = FALSE, include.deviance = FALSE,
                           include.groups = FALSE)

make_texreg("exp2_balancetests", list(ml2_7_extracted[[1]], ml2_7_extracted[[2]],
                                      ml2_38_extracted[[1]], ml2_38_extracted[[2]]),
            custom.header = list("Florida House District 7" = 1:2, "Florida House District 38" = 3:4),
            custom.model.names = c("Gloating Villain", "Report Card",
                                   "Gloating Villain", "Report Card"),
            custom.coef.map = list("bc_elections_voted" = "Past Elections Voted",
                                   "bc_mailaddressuspstype_H" = "Address Type (Apartment)",
                                   "age" = "Age",
                                   "black" = "Black",
                                   "hispanic" = "Hispanic",
                                   "asian" = "Asian",
                                   "female" = "Female",
                                   "married" = "Married",
                                   "bc_hh_size" = "Household Size",
                                   "(Intercept)" = "Constant"),
            custom.note = paste("\\item[\\hspace{-5mm}] %stars.",
                                "\\item[\\hspace{-5mm}] \\textit{Note:} Models estimated using multinomial logistic regression.
                                The dependent variable is each treatment condition regressed on the demographic covariate variables
                                used to balance the randomizations. The p-value of the likelihood ratio test was 0.76 for HD 7 and 
                                and 0.05 for HD 38, which means we fail to find that the covariates 
                                jointly predict treatment assignment."),
            multiple.tasks = "glm",
            scalebox = .8,
            caption = "Balance Tests of Treatment Conditions, Experiment 2 (FL)")







# *Experiments A and B -----------------------------------------------------------


# Tables 5 and 7. Means (pooled) ----------------------------------------------------------


emotions <- c("angry", "shame", "guilty", "disapp","indiff", "smug",  "happy", "proud")

means_pooled_expA <- map(emotions, ~ lm_robust(value ~ -1 + treatment, data = dfA,
                                           subset = emotion == .x & index != "net"))



make_texreg("expA_pooled", means_pooled_expA,
            stars = numeric(0),
            custom.header = list("Negative" = 1:4, 
                                 "Neutral" = 5,
                                 "Positive" = 6:8),
            custom.model.names = c("Angry", "Ashamed", "Guilty", "Disappointed", "Indifferent", 
                                   "Defiant/Smug", "Happy", "Proud"),
            custom.coef.map = list("treatmentGloating_Villain" = "Gloating Villain",
                                   "treatmentFoiled_Villain" = "Foiled Villain",
                                   "treatmentDisappointed_Hero" = "Disappointed Hero",
                                   "treatmentHappy_Hero" = "Happy Hero"),
            custom.note = paste("\\item[\\hspace{-5mm}] \\textit{Note:} The dependent variable is the emotion level on a 7-point Likert scale. 
            Emotions are ordered from negative, neutral, then positive. 
            In each column, the treatment that induces the highest level of emotion is bolded."),
            caption = "Mean Level of Anticipated Emotions by Treatment, Experiment A",
            resize.width = TRUE,
            refresh = FALSE,
            multiple.tasks = "glm")

map_dfr(1:8, ~means_pooled_expA[[.x]]$coefficients) %>% 
  select(treatmentGloating_Villain, treatmentFoiled_Villain,
         treatmentDisappointed_Hero, treatmentHappy_Hero) %>% 
  t() %>% 
  apply(2, function(x) x == max(x)) %>% 
  View()

means_netvote_expA <- map(emotions, ~ lm_robust(value ~ -1+ treatment, data = dfA,
                                            subset = emotion == .x & index == "net"))

make_texreg("expA_netvote", means_netvote_expA,
            custom.model.names = c("Angry", "Ashamed", "Guilty", "Disappointed", "Indifferent", 
                                   "Defiant/Smug", "Happy", "Proud"),
            custom.header = list("Negative" = 1:4, 
                                 "Neutral" = 5,
                                 "Positive" = 6:8),
            stars = numeric(0),
            custom.coef.map = list("treatmentGloating_Villain" = "Gloating Villain",
                                   "treatmentFoiled_Villain" = "Foiled Villain",
                                   "treatmentDisappointed_Hero" = "Disappointed Hero",
                                   "treatmentHappy_Hero" = "Happy Hero"),
            custom.note = paste("\\item[\\hspace{-5mm}] \\textit{Note:} The dependent variable is the difference in emotion levels, which ranges from -6 to 6. 
            Emotions are ordered from negative, neutral, then positive. 
            For negative and neutral (positive) emotions, the treatment that decreases (increases) the emotion level the most when voting is bolded."),
            caption = "Mean Effect of Voting Minus Not Voting on Level of Anticipated Emotion Levels by Treatment, Experiment A",
            multiple.tasks = "glm",
            refresh = FALSE,
            resize.width = TRUE)

map_dfr(1:8, ~means_netvote_expA[[.x]]$coefficients) %>% 
  select(treatmentGloating_Villain, treatmentFoiled_Villain,
         treatmentDisappointed_Hero, treatmentHappy_Hero) %>% 
  t() %>% 
  apply(2, function(x) x == min(x)) %>% 
  View()

## Note: when manually doing this, it's actually max for the positive emotions.


# Tables 6 and 8. Means (differences) ----------------------------------------------------------


emotions <- c("angry", "shame", "guilty", "smug", "happy", "proud")

means_pooled_expB <- map(emotions, ~ lm_robust(value ~ -1 + treatment, data = dfB,
                                           subset = emotion == .x & index != "net"))

make_texreg("expB_pooled", means_pooled_expB,
            custom.header = list("Negative" = 1:3, 
                                 "Positive" = 4:6),
            custom.model.names = c("Angry", "Ashamed", "Guilty",
                                   "Defiant/Smug", "Happy", "Proud"),
            custom.coef.map = list("treatmentGloating_Villain" = "Gloating Villain",
                                   "treatmentFoiled_Villain" = "Foiled Villain",
                                   "treatmentDisappointed_Hero" = "Disappointed Hero",
                                   "treatmentHappy_Hero" = "Happy Hero"),  
            stars = numeric(0),
            custom.note = paste("\\item[\\hspace{-5mm}] \\textit{Note:} 
                                The dependent variable is the emotion level on a 7-point Likert scale. 
            In each column, the treatment that induces the highest level of emotion is bolded."),      
            caption = "Mean Level of Anticipated Emotions by Treatment, Experiment B",
            multiple.tasks = "glm",
            refresh = FALSE,
            resize.width = TRUE)

map_dfr(1:6, ~means_pooled_expB[[.x]]$coefficients) %>% 
  select(treatmentGloating_Villain, treatmentFoiled_Villain,
         treatmentDisappointed_Hero, treatmentHappy_Hero) %>% 
  t() %>% 
  apply(2, function(x) x == max(x)) %>% 
  View()

means_netvote_expB <- map(emotions, ~ lm_robust(value ~ -1+ treatment, data = dfB,
                                            subset = emotion == .x & index == "net"))

make_texreg("expB_netvote", means_netvote_expB,
            stars = numeric(0),
            custom.header = list("Negative" = 1:3, 
                                 "Positive" = 4:6),
            custom.model.names = c("Angry", "Ashamed", "Guilty",
                                   "Defiant/Smug", "Happy", "Proud"),
            custom.coef.map = list("treatmentGloating_Villain" = "Gloating Villain",
                                   "treatmentFoiled_Villain" = "Foiled Villain",
                                   "treatmentDisappointed_Hero" = "Disappointed Hero",
                                   "treatmentHappy_Hero" = "Happy Hero"),
            custom.note = paste("\\item[\\hspace{-5mm}] \\textit{Note:} The dependent variable is the difference in emotion levels, which ranges from -6 to 6. 
            For negative (positive) emotions, the treatment that decreases (increases) the emotion level the most when voting is bolded."),
            caption = "Mean Effect of Voting Minus Not Voting on Level of Anticipated Emotions by Treatment, Experiment B",
            refresh = FALSE,
            multiple.tasks = "glm",
            resize.width = TRUE)

map_dfr(1:6, ~means_netvote_expB[[.x]]$coefficients) %>% 
  select(treatmentGloating_Villain, treatmentFoiled_Villain,
         treatmentDisappointed_Hero, treatmentHappy_Hero) %>% 
  t() %>% 
  apply(2, function(x) x == max(x)) %>% 
  View()

## Same as above, this time the maxes are harder to see so we did max instead of min



# Tables C1 and C2. Sample demographics ----------------------------------------------------------------
expA_sample_df <- read_csv("../original_data/Villain%2FHero+2+%2B+Policy+Ranking+Pilot+-+MTurk_June+6%2C+2017_14.10.csv") |> 
  slice(4:n()) |> 
  select(birthyear_1, contains("race"), gender, educ,
         contains("party")) |> 
  mutate(age = 2017-as.numeric(birthyear_1),
         female = case_when(gender == "Female" ~ 1,
                            TRUE ~ 0),
         education = case_when(educ == "Less than high school" ~ 1,
                               educ == "High school graduate, GED, or equivalent" ~ 2,
                               educ == "Some college" ~ 3,
                               educ == "2 year college degree" ~ 4,
                               educ == "4 year college degree" ~ 5,
                               educ == "Post-graduate degree" ~ 6),
         pid7 = case_when(party_str_rep == "Strong Republican" ~ 1,
                          party_str_rep == "Not a very strong Republican" ~ 2,
                          party_ind == "Republican" ~ 3,
                          party_ind == "No preference" | party_ind == "Don't know" ~ 4,
                          party_ind == "Democrat" ~ 5,
                          party_str_dem == "Not a very strong Democrat" ~ 6,
                          party_str_dem == "Strong Democrat" ~ 7)) |> 
  unite(race, contains("race"), na.rm = TRUE) |> 
  mutate(white = case_when(race == "White" ~ 1,
                           TRUE ~ 0)) |> 
  select(age:white)

expA_sample <- as.data.frame(round(apply(expA_sample_df, 2, mean, na.rm = TRUE), digits = 3)) %>% 
  cbind(round(apply(expA_sample_df, 2, sd, na.rm = TRUE), digits = 3)) %>% 
  cbind(round(apply(expA_sample_df, 2, min, na.rm = TRUE), digits = 3)) %>% 
  cbind(round(apply(expA_sample_df, 2, max, na.rm = TRUE), digits = 3)) %>% 
  rownames_to_column() %>% 
  mutate(rowname = case_when(rowname == "age" ~ "Age",
                             rowname == "female" ~ "Female",
                             rowname == "education" ~ "Education",
                             rowname == "pid7" ~ "Partisanship (7-point)",
                             rowname == "white" ~ "White"))
                             

make_kable(expA_sample, col.names = c("Variable", "Mean", "SD", "Min", "Max"),
           caption = "Summary Statistics, Experiment A")

expB_sample_df <- read_csv("../original_data/Govt_Spending__SSI__FINAL.csv") %>% 
  slice(2:n()) |> 
  filter(vh_js == 1) |> 
  select(contains("ps")) |> 
  rename(indinc = `ps2...53`,
         hhinc = `ps3...54`,
         education = ps0,
         race = `ps2...225`,
         birthyr = `ps3...228`,
         gender = ps4,
         party = ps5,
         party_ind = ps6,
         ideology = ps7) |> 
  mutate(female = case_when(gender == 2 ~1,
                            gender == 1 ~ 0),
         white = case_when(is.na(race) ~ NA_real_,
                           race == 1 ~ 1,
                           TRUE ~ 0),
         pid5 = case_when(party == 2 ~ 1, #Republican
                          party_ind == 2 ~ 2,
                          party_ind == 3 ~ 3,
                          party_ind == 1 ~ 4,
                          party == 1 ~ 5), #Democrat
         age = 2017-(as.numeric(birthyr)+1919)) |>  
  select(-contains("ps"), -c(race, party, party_ind, birthyr, gender)) |> 
  mutate(across(everything(), ~ as.numeric(.)))

expB_sample <- as.data.frame(round(apply(expB_sample_df, 2, mean, na.rm = TRUE), digits = 3)) %>% 
  cbind(round(apply(expB_sample_df, 2, sd, na.rm = TRUE), digits = 3)) %>% 
  cbind(round(apply(expB_sample_df, 2, min, na.rm = TRUE), digits = 3)) %>% 
  cbind(round(apply(expB_sample_df, 2, max, na.rm = TRUE), digits = 3)) %>% 
  rownames_to_column() %>% 
  mutate(rowname = case_when(rowname == "indinc" ~ "Individual Income",
                             rowname == "hhinc" ~ "Household Income",
                             rowname == "education" ~ "Education",
                             rowname == "ideology" ~ "Ideology",
                             rowname == "female" ~ "Female",
                             rowname == "pid5" ~ "Partisanship (5-point)",
                             rowname == "white" ~ "White",
                             rowname == "age" ~ "Age"))


make_kable(expB_sample, col.names = c("Variable", "Mean", "SD", "Min", "Max"),
           caption = "Summary Statistics, Experiment B")


# Tables C3 and C4 Means condition on voting --------------------------------------------------------
emotions <- c("angry", "shame", "guilty", "disapp","indiff", "smug",  "happy", "proud")

means_voting_expA <- map(emotions, ~ lm_robust(value ~-1 + treatment, data = dfA,
                                           subset = emotion == .x & index != "net" &
                                          vote == 1))


make_texreg("expA_voting", means_voting_expA,
            stars = numeric(0),
            custom.header = list("Negative" = 1:4, 
                                 "Neutral" = 5,
                                 "Positive" = 6:8),
            custom.model.names = c("Angry", "Ashamed", "Guilty", "Disappointed", "Indifferent", 
                                   "Defiant/Smug", "Happy", "Proud"),
            custom.coef.map = list("treatmentGloating_Villain" = "Gloating Villain",
                                   "treatmentFoiled_Villain" = "Foiled Villain",
                                   "treatmentDisappointed_Hero" = "Disappointed Hero",
                                   "treatmentHappy_Hero" = "Happy Hero"),
            custom.note = paste("\\item[\\hspace{-5mm}] \\textit{Note:} The dependent variable is the emotion level on a 7-point Likert scale. 
            Emotions are ordered from negative, neutral, then positive. 
            In each column, the treatment that induces the highest level of emotion is bolded."),
            caption = "Mean Level of Anticipated Emotions by Treatment When Voting, Experiment A",
            resize.width = TRUE,
            refresh = FALSE,
            multiple.tasks = "glm")

map_dfr(1:8, ~means_voting_expA[[.x]]$coefficients) %>% 
  select(treatmentGloating_Villain, treatmentFoiled_Villain,
         treatmentDisappointed_Hero, treatmentHappy_Hero) %>% 
  t() %>% 
  apply(2, function(x) x == max(x)) %>% 
  View()

emotions <- c("angry", "shame", "guilty", "smug", "happy", "proud")

means_voting_expB <- map(emotions, ~ lm_robust(value ~ -1 + treatment, data = dfB,
                                           subset = emotion == .x & index != "net" &
                                             vote == 1))

make_texreg("expB_voting", means_voting_expB,
            custom.header = list("Negative" = 1:3, 
                                 "Positive" = 4:6),
            custom.model.names = c("Angry", "Ashamed", "Guilty",
                                   "Defiant/Smug", "Happy", "Proud"),
            custom.coef.map = list("treatmentGloating_Villain" = "Gloating Villain",
                                   "treatmentFoiled_Villain" = "Foiled Villain",
                                   "treatmentDisappointed_Hero" = "Disappointed Hero",
                                   "treatmentHappy_Hero" = "Happy Hero"),  
            stars = numeric(0),
            custom.note = paste("\\item[\\hspace{-5mm}] \\textit{Note:} 
                                The dependent variable is the emotion level on a 7-point Likert scale. 
            In each column, the treatment that induces the highest level of emotion is bolded."),      
            caption = "Mean Levels of Anticipated Emotions by Treatment when Voting, Experiment B",
            multiple.tasks = "glm",
            refresh = TRUE,
            resize.width = TRUE)

map_dfr(1:6, ~means_voting_expB[[.x]]$coefficients) %>% 
  select(treatmentGloating_Villain, treatmentFoiled_Villain,
         treatmentDisappointed_Hero, treatmentHappy_Hero) %>% 
  t() %>% 
  apply(2, function(x) x == max(x)) %>% 
  View()


# Analysis for Experiments A and B ----------------------------------------
## Difference between Anger in Gloating Villain and Foiled Villain, Experiment A
summary(glht(means_pooled_expA[[1]], linfct= c("treatmentGloating_Villain - treatmentFoiled_Villain = 0")))

## Difference between Anger in Gloating Villain and Foiled Villain, Experiment B
summary(glht(means_pooled_expB[[1]], linfct = c("treatmentGloating_Villain - treatmentFoiled_Villain = 0")))

## Difference between Guilt in Gloating Villain and Disappointed Hero, Experiment B
summary(glht(means_pooled_expB[[3]], linfct = c("treatmentGloating_Villain - treatmentDisappointed_Hero = 0")))

## Net Difference between Anger in Gloating Villain and Happy Hero, Experiment A
summary(glht(means_netvote_expA[[1]], linfct = c("treatmentGloating_Villain - treatmentHappy_Hero = 0")))

## Net Difference between Anger in Gloating Villain and Happy Hero, Experiment A
summary(glht(means_netvote_expB[[1]], linfct = c("treatmentGloating_Villain - treatmentHappy_Hero = 0")))






# * Experiment C ----------------------------------------------------------


# Table D1. Sample demographics -------------------------------------------

expC_sample <- as.data.frame(round(apply(dfC_respondents |> select(-respondent_id), 2, mean, na.rm = TRUE), digits = 3)) %>% 
  cbind(round(apply(dfC_respondents |> select(-respondent_id), 2, sd, na.rm = TRUE), digits = 3)) %>% 
  cbind(round(apply(dfC_respondents |> select(-respondent_id), 2, min, na.rm = TRUE), digits = 3)) %>% 
  cbind(round(apply(dfC_respondents |> select(-respondent_id), 2, max, na.rm = TRUE), digits = 3)) %>% 
  rownames_to_column() |> 
  mutate(rowname = case_when(rowname == "d_age" ~ "Age",
                             rowname == "d_educ" ~ "Education",
                             rowname == "d_iswhite" ~ "White",
                             rowname == "d_isblack" ~ "Black",
                             rowname == "d_isasian" ~ "Asian",
                             rowname == "d_isother" ~ "Other",
                             rowname == "d_hispanic" ~ "Hispanic",
                             rowname == "d_female" ~ "Female",
                             rowname == "d_nonbinary" ~ "Non-binary",
                             rowname == "d_income" ~ "Income (refused = 6)",
                             rowname == "d_incomerf" ~ "Income refused",
                             rowname == "d_ideology" ~ "Ideology (Liberal to Conservative)",
                             rowname == "d_pid7" ~ "Partisanship (Democrat to Republican)"))

make_kable(expC_sample, col.names = c("Variable", "Mean", "SD", "Min", "Max"),
           caption = "Summary Statistics, Experiment C")

# Table D2. Current Emotions -----------------------------------------------

emotions <- c("angry", "shame", "guilty", "smug", "happy", "proud")

## Current Emotions
means_current_expC <- map(emotions, ~ lm_robust(value ~ Treatment, data = dfC,
                                        subset = emotion == .x & type == "current"))

make_texreg("expC_current", means_current_expC,
            custom.header = list("Negative" = 1:3, 
                                 "Positive" = 4:6),
            custom.model.names = c("Angry", "Ashamed", "Guilty",
                                   "Defiant/Smug", "Happy", "Proud"),
            custom.coef.map = list("TreatmentGloatingVillain" = "Gloating Villain",
                                   "TreatmentWorkVillain" = "Work Villain",
                                   "TreatmentPoliticalAngerReflection" = "Political Anger Reflection",
                                   "TreatmentAngerReflection" = "Non-Political Anger Reflection",
                                   "(Intercept)" = "Constant (Baseline = Control)"),  
            custom.note = paste("\\item[\\hspace{-5mm}] %stars.",
                                "\\item[\\hspace{-5mm}] \\textit{Note:} 
                                Models estimated using ordinary least squares regression,
                                with robust standard errors.
                                The dependent variable is the current emotion level on a 7-point Likert scale."),      
            caption = "Effect of Treatments on Levels of Current Emotions, Experiment C",
            multiple.tasks = "glm",
            refresh = TRUE,
            resize.width = TRUE)




# Table D3. Voting Emotions ------------------------------------------------

## Anticipated Emotions when Voting
means_vote_expC <- map(emotions, ~ lm_robust(value ~ Treatment, data = dfC,
                                                subset = emotion == .x & type == "vote"))

make_texreg("expC_vote", means_vote_expC,
            custom.header = list("Negative" = 1:3, 
                                 "Positive" = 4:6),
            custom.model.names = c("Angry", "Ashamed", "Guilty",
                                   "Defiant/Smug", "Happy", "Proud"),
            custom.coef.map = list("TreatmentGloatingVillain" = "Gloating Villain",
                                   "TreatmentWorkVillain" = "Work Villain",
                                   "TreatmentPoliticalAngerReflection" = "Political Anger Reflection",
                                   "TreatmentAngerReflection" = "Non-Political Anger Reflection",
                                   "(Intercept)" = "Constant (Baseline = Control)"),  
            custom.note = paste("\\item[\\hspace{-5mm}] %stars.",
                                "\\item[\\hspace{-5mm}] \\textit{Note:} 
                                Models estimated using ordinary least squares regression,
                                with robust standard errors.
                                The dependent variable is the anticipated emotion level when voting on a 7-point Likert scale."),      
            caption = "Effect of Treatments on Levels of Anticipated Emotions When Voting, Experiment C",
            multiple.tasks = "glm",
            refresh = TRUE,
            resize.width = TRUE)


# Table D4. Not Voting Emotions ------------------------------------------------


## Anticipated Emotions when Not Voting
means_novote_expC <- map(emotions, ~ lm_robust(value ~ Treatment, data = dfC,
                                             subset = emotion == .x & type == "novote"))

make_texreg("expC_novote", means_novote_expC,
            custom.header = list("Negative" = 1:3, 
                                 "Positive" = 4:6),
            custom.model.names = c("Angry", "Ashamed", "Guilty",
                                   "Defiant/Smug", "Happy", "Proud"),
            custom.coef.map = list("TreatmentGloatingVillain" = "Gloating Villain",
                                   "TreatmentWorkVillain" = "Work Villain",
                                   "TreatmentPoliticalAngerReflection" = "Political Anger Reflection",
                                   "TreatmentAngerReflection" = "Non-Political Anger Reflection",
                                   "(Intercept)" = "Constant (Baseline = Control)"),  
            custom.note = paste("\\item[\\hspace{-5mm}] %stars.",
                                "\\item[\\hspace{-5mm}] \\textit{Note:} 
                                Models estimated using ordinary least squares regression,
                                with robust standard errors.
                                The dependent variable is the anticipated emotion level when not voting on a 7-point Likert scale."),      
            caption = "Effect of Treatments on Levels of Anticipated Emotions When Not Voting, Experiment C",
            multiple.tasks = "glm",
            refresh = TRUE,
            resize.width = TRUE)


# Analysis for Experiment C -----------------------------------------------

# 1st paragraph: When voting, comparing gloating villain with work villain treatments

## Anger, GV
summary(glht(means_vote_expC[[1]], linfct = c("TreatmentGloatingVillain = 0")))

### one-tailed
summary(glht(means_vote_expC[[1]], linfct = c("TreatmentGloatingVillain = 0")))$test$pvalues/2

## Anger, GV - WV
summary(glht(means_vote_expC[[1]], linfct = c("TreatmentGloatingVillain - TreatmentWorkVillain = 0")))

### one-tailed
summary(glht(means_vote_expC[[1]], linfct = c("TreatmentGloatingVillain - TreatmentWorkVillain = 0")))$test$pvalues/2

## Smug, GV
summary(glht(means_vote_expC[[4]], linfct = c("TreatmentGloatingVillain = 0")))

### one-tailed
summary(glht(means_vote_expC[[4]], linfct = c("TreatmentGloatingVillain = 0")))$test$pvalues/2

## Smug, GV - WV
summary(glht(means_vote_expC[[4]], linfct = c("TreatmentGloatingVillain - TreatmentWorkVillain = 0")))

### one-tailed
summary(glht(means_vote_expC[[4]], linfct = c("TreatmentGloatingVillain - TreatmentWorkVillain = 0")))$test$pvalues/2

# 2nd paragraph: Effect of GV on current emotional states

## Anger, GV
summary(glht(means_current_expC[[1]], linfct = c("TreatmentGloatingVillain = 0")))

###one-tailed
summary(glht(means_current_expC[[1]], linfct = c("TreatmentGloatingVillain = 0")))$test$pvalues/2

## Smugness, GV
summary(glht(means_current_expC[[4]], linfct = c("TreatmentGloatingVillain = 0")))

### one-tailed
summary(glht(means_current_expC[[4]], linfct = c("TreatmentGloatingVillain = 0")))$test$pvalues/2


# 3rd paragraph: Next smallest p-values
## It's guilt for Current
summary(glht(means_current_expC[[3]], linfct = c("TreatmentGloatingVillain = 0")))


### one-tailed
summary(glht(means_current_expC[[3]], linfct = c("TreatmentGloatingVillain = 0")))$test$pvalues/2

## And it's shame for Voting
summary(glht(means_vote_expC[[2]], linfct = c("TreatmentGloatingVillain = 0")))

### one-tailed
summary(glht(means_vote_expC[[2]], linfct = c("TreatmentGloatingVillain = 0")))$test$pvalues/2


## But shame is not very high for not voting
summary(glht(means_novote_expC[[2]], linfct = c("TreatmentGloatingVillain = 0")))

### one-tailed
summary(glht(means_novote_expC[[2]], linfct = c("TreatmentGloatingVillain = 0")))$test$pvalues/2

## and also guilty for Voting
summary(glht(means_vote_expC[[3]], linfct = c("TreatmentGloatingVillain = 0")))

### one-tailed
summary(glht(means_vote_expC[[3]], linfct = c("TreatmentGloatingVillain = 0")))$test$pvalues/2

# 4th paragraph: Finally, we compare GV to other treatments.

## GV vs. Non-Political Anger Reflection for Smug
summary(glht(means_vote_expC[[4]], linfct = c("TreatmentGloatingVillain - TreatmentAngerReflection = 0")))

### one-tailed
summary(glht(means_vote_expC[[4]], linfct = c("TreatmentGloatingVillain - TreatmentAngerReflection = 0")))$test$pvalues/2

## GV vs. Political Anger Reflection for Smug
summary(glht(means_vote_expC[[4]], linfct = c("TreatmentGloatingVillain - TreatmentPoliticalAngerReflection = 0")))

### one-tailed
summary(glht(means_vote_expC[[4]], linfct = c("TreatmentGloatingVillain - TreatmentPoliticalAngerReflection = 0")))$test$pvalues/2

# Misc reviewer memo analysis

## Difference between GV and WV for current anger
summary(glht(means_current_expC[[1]], linfct = c("TreatmentGloatingVillain - TreatmentWorkVillain = 0")))


### one-tailed
summary(glht(means_current_expC[[1]], linfct = c("TreatmentGloatingVillain - TreatmentWorkVillain = 0")))$test$pvalues/2

## Difference between GV and WV for current smug
summary(glht(means_current_expC[[4]], linfct = c("TreatmentGloatingVillain - TreatmentWorkVillain = 0")))

## Difference between PAR and GV not significant for anticipated anger
summary(glht(means_vote_expC[[1]], linfct = c("TreatmentGloatingVillain - TreatmentPoliticalAngerReflection = 0")))


# Table E1 ----------------------------------------------------------------
heroes <- read_csv("Table_Heroes_By_Party.csv") |> 
  # Create labels for each row for Dem, Ind, or Rep
  mutate(subject_party = case_when(str_detect(Hero, "Democrat") ~ "Democrat",
                                   str_detect(Hero, "Republican")~ "Reublican",
                                   str_detect(Hero, "Independent") ~ "Independent")) |> 
  # Create a total_n so that I can create percentages column
  mutate(total_n = case_when(str_detect(Hero, "Democrat") ~ str_extract(Hero, "\\d\\d\\d"),
                             str_detect(Hero, "Republican")~ str_extract(Hero, "\\d\\d\\d"),
                             str_detect(Hero, "Independent") ~ str_extract(Hero, "\\d\\d\\d"))) |> 
  # Fill down
  fill(subject_party, total_n) |> 
  # Create percentages
  mutate(`Percent of Total` = paste0(round(N/as.numeric(total_n)*100), "\\%")) |> 
  mutate(`Percent of Total` = case_when(is.na(Party) ~ NA_character_,
                                        TRUE ~ `Percent of Total`)) |> 
  # Pick top 10 for each
  group_by(subject_party) |> 
  slice(1:11) |> 
  ungroup() |> 
  select(-c(subject_party, total_n))

make_kable(heroes, col.names = names(heroes),
           resize = TRUE,
           caption = "List of Heroes by Subject's Party Identification")

villains <- read_csv("Table_Villains_By_Party.csv") |> 
  # Create labels for each row for Dem, Ind, or Rep
  mutate(subject_party = case_when(str_detect(Villain, "Democrat") ~ "Democrat",
                                   str_detect(Villain, "Republican")~ "Reublican",
                                   str_detect(Villain, "Independent") ~ "Independent")) |> 
  # Create a total_n so that I can create percentages column
  mutate(total_n = case_when(str_detect(Villain, "Democrat") ~ str_extract(Villain, "\\d\\d\\d"),
                             str_detect(Villain, "Republican")~ str_extract(Villain, "\\d\\d\\d"),
                             str_detect(Villain, "Independent") ~ str_extract(Villain, "\\d\\d\\d"))) |> 
  # Fill down
  fill(subject_party, total_n) |> 
  # Create percentages
  mutate(`Percent of Total` = paste0(round(N/as.numeric(total_n)*100), "\\%")) |> 
  mutate(`Percent of Total` = case_when(is.na(Party) ~ NA_character_,
                                        TRUE ~ `Percent of Total`)) |> 
  # Pick top 10 for each
  group_by(subject_party) |> 
  slice(1:11) |> 
  ungroup() |> 
  select(-c(subject_party, total_n))

make_kable(villains, col.names = names(villains),
           resize = TRUE,
           caption = "List of Villains by Subject's Party Identification")
         
